/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of single_lambda_emergence methods.

// $Id$

#include "single_lambda_emergence.h"
#include "emergence-fits.h"
#include "CCfits/CCfits"
#include "sstream"
#include "biniostream.h"
#include "fits-utilities.h"
#include "constants.h"

using namespace CCfits;
using namespace blitz;

/** Constructor sets up the emergence object from a FITS file.  All
    functionality is in the emergence class.  */
mcrx::single_lambda_emergence::single_lambda_emergence (const CCfits::FITS& file):
    emergence<T_float> (file)
{};

/** Creates FITS data cubes for saving the camera image arrays.  The
    number of slices (i.e., wavelengths) is specified. \bug Doesn't
    work for 1 slice. */
void mcrx::single_lambda_emergence::create_HDU (FITS& file, int slices,
						const T_unit_map& units) const
{
  // convert image values (which is flux in the pixels with length
  // units the same as the camera distance units) to surface
  // brightness. This involves converting the 1/area unit to 1/m and
  // dividing by the solid angle of the pixels. We use units to
  // convert whatever length unit we have to m.
  const T_float to_m = units::convert(units.get("length"),"m");
  const string sb_unit = units.get("L_lambda") + "/m^2/sr";
  const string sb_factr_unit = "m^2/(" + units.get("length") + "^2*sr)";
  // (1kpc^2 = 9.5214e38 m^2)

  for (T_camera_map::const_iterator i = camera_numbers.begin();
       i != camera_numbers.end(); ++i) {

    const int camera_number = i->first;
    const T_camera& c = get_camera (camera_number);

    const array_2& image = get_image (camera_number);
    assert (image.size()> 0);
    std::ostringstream ost;
    ost << "CAMERA" << camera_number << suffix (); 

    // We only want to create the HDU's if they don't already exist,
    ExtHDU* hdu;
    try {
      hdu = &open_HDU (file,ost.str());
    }
    catch (CCfits::FITS::NoSuchHDU&) {
      std::vector<long> naxes;
      naxes.push_back(image.extent(firstDim ) );
      naxes.push_back(image.extent(secondDim ) );
      naxes.push_back(slices);
      hdu = file.addImage (ost.str(), FLOAT_IMG, naxes);
      //std::cout << "creating HDU " << ost.str()<< endl;

      hdu->writeComment(comment ());

      // write wcs keywords
      c.write_wcs(*hdu, units);

      const T_float sb_factr = 1./(to_m*to_m*c.pixel_solid_angle());
      assert(sb_factr > 0);

      // write units
      hdu->addKey("imunit", sb_unit,
		  "Pixel quantity is surface brightness");
      hdu->addKey("pixel_sr", c.pixel_solid_angle(),
		  "[sr] Solid angle of pixels");
      hdu->addKey("sb_factr", sb_factr,"["+sb_factr_unit +
		  "] Conversion from internal camera flux to surface brightness");
    }
  } 
  
}

/** Writes the images to a slice in the data cubes in HDU's
    CAMERAi-suffix ().  The HDU's must have been created with a
    previous call to write_HDU. For the file format specification see
    the mcrx-readme. */
void mcrx::single_lambda_emergence::write_images (CCfits::FITS& file,
						  T_float normalization,
						  int slice) const
{
  assert(normalization > 0);

  cout << "Writing " << n_cameras_loaded() << " cameras: \n";

  for (T_camera_map::const_iterator i = camera_numbers.begin();
       i != camera_numbers.end(); ++i) {
    const int camera_number = i->first;

    // find HDU
    std::ostringstream ost;
    ost << "CAMERA" << camera_number << suffix ();
    cout << ost.str() << endl;

    ExtHDU& hdu = open_HDU (file, ost.str ());

    std::ostringstream ost2;
    ost2 << "normalization" << slice;
    hdu.addKey(ost2.str(), normalization,
	       "Image normalization (internal usage)");

    double sb_factr = 1; // for some reason CCfits makes a double keyword!?
    try {
      hdu.readKey ("sb_factr", sb_factr);
    }
    catch (CCfits::HDU::NoSuchKeyword&) {}

    // Write image array 
    // We have storage order issues here, FITS images are column major
    const array_2& image = get_image (camera_number);
    assert (all (image < blitz::huge (T_float ())));
    assert (all (image >= 0));
    array_2 temp_image (image.shape(), ColumnMajorArray<2> ());
    assert(temp_image.isStorageContiguous());

    temp_image = (normalization*sb_factr)*image;
    std::valarray<T_float> temp (temp_image.dataFirst(), temp_image.size());

    // write the proper slice of the data cube
    std::vector<long> begin, end;
    begin.push_back(1 );
    begin.push_back(1 );
    begin.push_back(slice + 1);
    end.push_back(temp_image.extent(firstDim ) );
    end.push_back(temp_image.extent(secondDim) );
    end.push_back(slice + 1  );
    int s = 1;
    for (vector<long>::const_iterator ii = begin.begin(), jj = end.begin(); 
	 ii != begin.end(); ++ii, ++jj) {
      s*= (*jj - *ii + 1);
    }
    assert (s == temp.size());

    hdu.write(begin, end, temp );
  }
}

/** Loads the camera image arrays from a previously saved FITS file. */    
void mcrx::single_lambda_emergence::load_images (const CCfits::FITS& file,
						 int slice)
{
  cout << "Loading " << n_cameras_loaded() << " cameras: \n";
  for (T_camera_map::iterator i = camera_numbers.begin();
       i != camera_numbers.end(); ++i) {

    const int camera_number = i->first;
    T_camera& c = get_camera (camera_number);

    std::ostringstream ost;
    ost << "CAMERA" << camera_number << suffix ();
    cout << ost.str() << endl;
    //const_cast<CCfits::FITS&> (file).read(string (ost.str()));
    const CCfits::ExtHDU& hdu = open_HDU (file,ost.str() );
    assert (hdu.axis(0) == c.xsize ());
    assert (hdu.axis(1 ) == c.ysize ());
    array_2& image = get_image (camera_number);

    //read the proper slice of the data cube
    std::vector<long> begin, end;
    begin.push_back(1 );
    begin.push_back(1 );
    begin.push_back(slice + 1);
    end.push_back(image.extent(firstDim ) );
    end.push_back(image.extent(secondDim) );
    end.push_back(slice + 1  );
    int s = 1;
    for (vector<long>::const_iterator ii = begin.begin(), jj = end.begin(); 
	 ii != begin.end(); ++ii, ++jj) {
      s*= (*jj -*ii + 1);
    }

    read(hdu, image);

    // check for a normalization keyword
    float normalization = 1 ;
    double sb_factr = 1;
    try {
      const_cast<CCfits::ExtHDU&> (hdu).readKey (string ("normalization"),
						 normalization);
    }
    catch (...) {}// CCfits::HDU::NoSuchKeyword&) {}
    try {
      const_cast<CCfits::ExtHDU&> (hdu).readKey (string ("sb_factr"),
						sb_factr);
    }
    catch (...) {}// CCfits::HDU::NoSuchKeyword&) {}

    image *= 1/(normalization*sb_factr);
  }
}


/** Dump the image arrays to a file.  The main objective here is to be
    fast, since we have gotten a termination signal. */
void mcrx::single_lambda_emergence::write_dump (binofstream& file) const
{
  for (T_camera_map::const_iterator i = camera_numbers.begin();
       i != camera_numbers.end(); ++i) {
    const int camera_number = i->first;

    const array_2& image = get_image (camera_number);
    // We strictly don't need to save the image size here but it's a good check
    array_2::T_index extent = image.shape();
    file << extent [0]
	 << extent [1]; 

    assert(image.isStorageContiguous());
    file.write(reinterpret_cast<const char*> (image.dataFirst()),
	       image.storageSize()*sizeof (array_2::T_numtype) );
  } 
}

/** Read the image arrays from a previously written dump file. */
void mcrx::single_lambda_emergence::load_dump (binifstream& file)
{
  for (T_camera_map::iterator i = camera_numbers.begin();
       i != camera_numbers.end(); ++i) {

    const int camera_number = i->first;
    T_camera& c = get_camera (camera_number);

    array_2& image = get_image (camera_number);
    array_2::T_index size;
    file >> size [0]
	 >>  size [1]; 
    //cout << size << endl;
    assert (size [0] == c.xsize ());
    assert (size [1] == c.ysize ());

    assert(image.isStorageContiguous());
    file.read (reinterpret_cast<char*> (image.dataFirst()),
	       image.storageSize()*sizeof (array_2::T_numtype) );
  }
}


